A model for the large-scale circulation in turbulent Rayleigh-Benard convection 
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A model for the large-scale circulation (LSC) dynamics of turbulent Rayleigh-Benard convection 
is presented. It consists of two stochastic ordinary differential equations motivated by the Navier- 
Stokes equation, one each for the strength and the azimuthal orientation of the LSC. Stochastic 
forces represent phenomenologically the action of the turbulent fluctuations on the LSC. Without 
adjustable parameters, the model yields a meandering LSC with occasional rotations, and with 
more rare cessations that occur a few times per day - differing only by about a factor of two 
from experiment. Also as in experiments, the distribution of LSC orientation-change is uniform for 
cessations and a power law for rotations. 

PACS numbers: 47.27.-i, 05.65.+b, 47.27.te, 47.27.eb 



Rayleigh Benard convection (RBC) consists of a fluid- 
filled container heated from below It is well-known 
that cylindrical containers with aspect ratio Y = D/ L ~ 
1 (D is the sample diameter and L the height) have a 
large-scale circulation (LSC), also known as the "mean 
wind" @, i, i, i, i, E S i- The LSC consists of 
warm up-flowing and cool down-flowing fluid on opposite 
sides of the sample. Its near-vertical circulation plane 
has an orientation Oq that undergoes azimuthal diffusion 
Mini- a n interesting aspect is that 8q, in addition to 
its diffusive meandering, undergoes relatively rapid os- 
cillations 0] and, on somewhat longer time scales, expe- 
riences spontaneous and erratic re-orientations through 
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an azimuthal displacement AO. 
One mechanism for re-orientation is an azimuthal rota- 
tion of the entire structure without much change of the 
flow speed. [13] Rotations lead to a power-law distribu- 
tion of AO, with small rotations more likely than large 
ones. [ToL EH The other mechanism is a cessation of the 
flow in which it stops briefly and then starts up again 
in a random new orientation, resulting in a uniform dis- 
IcL T3 |. Aside from its fundamental in- 



tribution of AO 
terest, the phenomenon of reorientation is important for 
instance because it occurs in natural convection of the at- 



mosphere [16J, and because it is responsible for changes 
in the orientation of Earth's magnetic field when it occurs 
in Earth's outer core 



Several models were proposed to reproduce the LSC 
dynamics. Those of Sreenivasan et al. and Benzi 
[l8| were based on similar stochastic equations that were 
chosen so as to produce the desired reversal events. The 
equations yielded two opposite stable flow directions with 
transitions between them, in qualitative agreement with 
experiments by Sreenivasan et al. 7]. However, such 
local flow reversals usually are not cessations or large ro- 
tations; they are equivalent to crossings of the LSC orien- 
tation past a fixed angle and correspond mostly to small 
orientation changes, or "jitter" [loj ] . Another model, by 
Fontenele Araujo et al. [lSj, was based on a deterministic 
force balance between buoyancy and drag. It produced 



LSC reversals (AO — n), but did not consider azimuthal 
motion and thus did not reproduce the azimuthal dynam- 
ics seen in experiments. Finally, a deterministic model 
based on the Boussinesq equations with slip boundary 
conditions and for ellipsoidal sample geometries was de- 
veloped by Resagk et al. [23]. It neglected dissipation 
and diffusion, but added these effects phenomenologically 
after the model was derived. The result was a set of ordi- 
nary differential equations (ODEs) with several parame- 
ters that could be tuned to produce a LSC with various 
dynamical states, including oscillations and chaotic me- 
andering; but to our knowledge the existence and statis- 
tics of rotations and cessations were not explored. 

We present a model consisting of two coupled stochas- 
tic ODEs, one for the speed (or "strength") of the LSC 
and the other for the azimuthal LSC orientation. We 
retained the physically important terms of the Navier- 
Stokes (NS) equation and took volume averages. The 
resulting deterministic ODE for the LSC strength rep- 
resents the balance between buoyancy and drag forces 
and has two fixed points; one stable and the other, corre- 
sponding to cessations, unstable. It differs from the equa- 
tions of Refs. [3] and [13] in that the nonlinearity is of 
order 3/2 rather than of the more usual cubic order. The 
second ODE describes the azimuthal motion of the LSC 
which is suppressed by its angular momentum. The dy- 
namics of the model arises from the addition of stochas- 
tic forces that represent in a phenomcnological sense the 
action of the turbulent fluctuations that exist through- 
out the system interior. We determined some parameters 
of the model from independent measurements that did 
not involve reorientations per se, determined others from 
theory, and did not adjust any of them arbitrarily. An 
interesting physical aspect of the model is that the an- 
gular momentum of the LSC hinders reorientations when 
the flow is vigorous; when the flow becomes feeble, then 
the angular momentum is small and the stochastic forces 
are able to cause significant orientation change. We note 
that the model of Resagk et al. [13] suggests that the 
angular momentum of the LSC also plays an important 
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role in its oscillations. Our model does not produce LSC 
oscillations, but yields cessations that are only slightly 
more frequent than the experimental observation of one 
or two per day (lol . [l4j ] . In agreement with experiment 
gives a uniform probability distribution p(A9) 
for cessations and a power-law distribution for rotations. 

For the LSC strength we consider the velocity com- 
ponent U0, where </> is an angle that sweeps the plane 
of the LSC, and describes it without azimuthal motion. 
One expects the acceleration to be due to a balance be- 
tween buoyancy and drag forces. Thus we include in the 
NS equation for only the acceleration, buoyancy, and 
viscous drag terms, and neglect the nonlinear term 21] : 



= ga(T - T ) + vS7 2 u 4 



(1) 



Here a is the isobaric thermal expansion coefficient, g the 
acceleration of gravity, and v the kinematic viscosity. 

To obtain a model in the form of an ODE that de- 
scribes the flow with only a few variables, we take a global 
average over the field variables that retains the essential 
physics. To carry out the a verage, we consider the exper- 
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14 . |22| that the temperature 



imental observation 
of the LSC at the side wall at mid-height can be written 
as 



T = T a + 5cos(6 -' 



(2) 



where the temperature amplitude S represents the 
strength of the LSC, and where 9q is its azimuthal ori- 
entation. The buoyancy acts on the entire LSC and is 
proportional S. The profile is taken to be given by Eq. [21 
and proportional to the cylindrical radius r. The velocity 
is assumed to be linear in r and a step function in 9. Note 
that these assumptions about the geometry of the flow 
only affect the numerical prcfactors in the equations, and 
not the functional form. The drag is assumed to occur in 
the viscous boundary layers, so V 2 u^, « U/X 2 (U is the 
maximum speed near the side wall, A is the boundary- 
layer width), where A = (L/2) x R~V 2 (R e>i = UL/v 
is the instantaneous Reynolds number). The volume av- 
erage requires another factor of QX/L since the drag is 
mainly in the boundary layers. These approximations 
result in the volume-averaged equation 



(2/3)f7 = (2/37r).ga<S - YlvUR^jL 



(3) 



Next we make the assumption that the amplitude 5 is 
instantaneously proportional to the speed U , since both 
variables are measures of the LSC strength. To determine 
the proportionality, we find the steady-state solution 



of Eq. [3] Here R e is the normal steady-state Reynolds 
number. In the dynamical equation, we allow the drag to 
depend on the instantaneous value R e ,i, so the drag term 
instantaneously scales as [7 3 / 2 . We substitute Eq. [4] into 
Eq. [31 combine all parameters into two constants, and 
add a noise term that represents the turbulent fluctua- 
tions of the flow to get the Langevin equation 



• S S 3 / 2 
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fs(t) 



with 



8q = l8irATaRl /2 /R ; t s = L 2 j{\%vR\ 



(5) 



(6) 



Here the Rayleigh number is R = agATL 3 / kv with AT 
the applied temperature difference and k the thermal dif- 
fusivity, and the Prandtl number is a = v j k. 

In the absence of the noise term Eq. [5] has two fixed 
points, one unstable at S = and one stable when 5 = 5q. 
In the stochastic equation this feature reproduces some 
of the dominating behavior of the LSC; the LSC spends 
most of its time meandering near the stable fixed point 
at 5q, but occasionally it ceases when fluctuations drive 
it close to S = 0. 

To study the predictions of the model, we consider the 
example R = 1.1 x 10 10 and a = 4.4 for a sample with 
L = 24.76 cm [13]. Measurements yielded S = 0.25 
K and R e = 3700 [Hj]. From Eq. M one has t s = 85 
s and 8q = 0.10 K. For So theory and experiment are 
in order of magnitude agreement, which is as much as 
we can expect given the approximations made in the 
model derivation. For numerical calculations we adopt 
the experimental value because it presumably is more 
appropriate for the physical system. To gain information 
about the noise intensity, we examined the experimental 
mean-square amplitude-change ((dS) 2 ) ((...) represents 
a time average) over a time period dt as a function of 
dt. For time scales that were not too large we found 
({dS) 2 } = D s dt with D 5 = 3.5 x 10- 5 K 2 /s, suggesting 
a diffusive process. This method was used before [ll( to 
determine the diffusivity of 9q] but in the present case the 
diffusive scaling holds only over intermediate time scales 
because S is bounded. With this experimental input we 
make the noise in the model Brownian with diffusivity 
Ds, so fs(t) is Gaussian distributed with width yjDg/h 
where h is the time step in the simulation. 

The frequency of cessations is given approximately by 
the Arrhenius-Kramers result for diffusion over a poten- 
tial barrier AV 24|. Here V = — J SdS, integrated over 
the deterministic part of Eq.El and AV = V(Q)— V(Sq) = 
8$/(10t s ) ~ 7.3 ■ 10~ 5 K 2 /s. When AV > D s , the rate 
of cessations u> is given by 



(2/3n)ga6 = 12vUR l J 2 /L 2 



(4) 



LV 



ujo exp(—AV/Ds 



(7) 
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where lu = 1/(2V2ttt s ) ~ 1.3 ■ 1CT 3 s" 1 0. This yields 
about 14 cessations per day, an order of magnitude more 
than the experimental value of one or two per day. We 
attribute this difference to the fact that the condition 
AV » Ds is not really satisfied. 

The second Langevin equation describes the azimuthal 
motion. The main driving force is the turbulent noise. 
We estimated that angular momentum of the LSC in the 
<f> coordinate damps the relatively slow azimuthal motion 
much more than the viscous drag across the boundary 
layer near the side wall. The physical explanation is that 
a rotating body has some stability due to its angular 
momentum and requires a larger torque to rotate in an 
orthogonal direction than a non-rotating body. This phe- 
nomenon is represented by the transport term in the NS 
equation. Thus, neglecting the drag term (and Earth's 
Coriolis force; see [LJJ), we have 



ue + (u ■ V)ug = 



(8) 



A volume average gives (1/3)L9q = — (2/3)U8q . The 
angular momentum leads to an effective damping that is 
proportional to the wind strength U, which is important 
for understanding the azimuthal dynamics during cessa- 
tions. Again, we convert U to 5, combine the remaining 
parameters to get a new constant, and add a noise term 
representing turbulent fluctuations to get 



?0 = 



taSo 



(9) 



with Tg = L 2 /{2vR e ) . For R = 1.1 X 10 10 and a = 4.4 
the model predicts r g - ~ 13 s. The turbulent noise in this 
coordinate is also found to be Brownian, with diffusivity 
D e - = 2.5 x 1CP 5 rad 2 /s 3 . This diffusivity comes from 
a fit of ((d9 ) 2 ) = Dgdt to experimental data for 9o{t) 
(((d9o) 2 } is the mean-square change in rotation rate over 
the time period dt). Again, this scaling only holds for 
intermediate time periods because 9q is bounded. 

The two stochastic ODEs Eqs. [5] and [9] are our model 
for the LSC dynamics. Using the experimentally deter- 
mined values of <5o, Ds and D g - and predictions for t$ and 
Tg based on the measured value of R e discussed above, 
they can be integrated to get time series for 5 and 9 . We 
used a simple first-order Euler method to solve the equa- 
tions stepwise with a time step shorter than the smallest 
timescale t - of the system. Figure Q] shows a simulated 
time series over about one day. One can see, as we ex- 
pected, that the LSC amplitude 5 is stable with an occa- 
sional cessation where the amplitude drops to zero. From 
much longer simulations we found that cessations occur 
about 3.8 times per day, which (we presume because the 
relation AV ^> Ds is not satisfied very well) is a factor 
of 3.7 less than the result from Eq.[7|and about twice the 
frequency measured experimentally [10]. Considering the 




FIG. 1: A time series of 5 and Qo from the simulation of 
Eqs. [S]and[5] The LSC strength S fluctuates around (5o, and 
occasionally those fluctuations are large enough to cause a 
cessation where 8 ~ 0. 




time (s) 



FIG. 2: A short section of the time series of S and 9o showing 
a cessation in detail. The amplitude 5 drops to near zero; 
while 5 is small, the azimuthal motion becomes fast because 
it is no longer suppressed by the LSC angular momentum. 



approximations made in the model derivation, we regard 
this as very satisfying agreement with experiment. 

Figure [T] shows that the orientation meanders as ex- 
pected, but one must look on a shorter time scale to 
see the details of the dynamics. Thus, Fig. [2] shows a 
shorter section of the same time series that contains a 
cessation. One sees how S gradually drops to zero, then 
grows back up again over a few hundred seconds, just as 
observed experimentally (hJ. The time series for 9q is 
interesting because there is a large change in 9q during 
the cessation, again as seen experimentally [lOj . Equa- 
tion [5] for the azimuthal motion implies that, when 5 and 
thus the angular momentum are small during cessations, 
the damping term becomes small so the turbulent fluctu- 
ations are free to accelerate the LSC to large azimuthal 
rotation rates. When the LSC is strong, i. e. 6 « 5q, then 
the larger angular momentum of the LSC supresses the 
azimuthal rotation. This inverse relationship between 6 
and 9q was observed in experiments 13, 14 1 but had not 
been explained by any previous model. 

In order to determine the statistics of cessations and 
rotations, we analyzed the simulated time scries using the 
same algorithms as those applied before to experimental 
data lCj, [lj] . Figure [3k. shows the probability distribu- 
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FIG. 3: (a)The probability distribution of the azimuthal 
change AO during cessations. The solid line represents the 
uniform distribution p(A8) — 1/n. (b) p(A8) for rotations. 
The solid line is a power law. 



tion p(A0) of the orientation change AO during cessa- 
tions. The results from the simulations are consistent 
with a uniform distribution p{A9) = in agreement 



with the experiments [10lll4| . This is an important result 
that no model had predicted before, either because the 
dependence was not contained in it, @, EH or because 
it was assumed that cessations would result in a reversal 
of the flow direction, i.e. AO = tt or because the 
issue was not addressed [ioj]. The angular momentum 
that usually suppresses the azimuthal motion of the LSC 
is reduced during cessations, allowing turbulent fluctu- 
ations to azimuthally rotate the LSC freely. For large 
enough noise strengths, the azimuthal distance traveled 
is large over the duration of the cessation, and this results 
in a final orientation independent of the orientation be- 
fore the cessation, which explains the uniform p(A0). We 
also find rotations in the model to result from a similar 
mechanism. Rotations typically occur when the angu- 
lar momentum of the LSC is still large, so the azimuthal 
rotation in limited. This results in a monotonically de- 
creasing distribution of AO for rotations. However, the 
azimuthal rotation rate can become large even when 5 
is somewhat lower than normal, resulting in more rota- 
tions, and in particular more large rotations, than would 
be expected from purely Brownian noise with constant 
damping coefficient The simulation results can be 

represented reasonably well by a power law for p(A0) as 
shown in Fig [3Jd, in qualitative agreement with the ex- 
periments [14| . 

To summarize, we derived a dynamical model moti- 
vated by the Navier-Stokes equations to describe the LSC 
in terms of two variables: its strength S and its azimuthal 
orientation Oq. The model used as input four parame- 
ters that were determined from independent experimen- 
tal measurements, and required no adjustable parame- 
ters. Each term in the equations has a clear physical 
meaning. The model produces a stable LSC with occa- 
sional cessations and rotations, and p{ AO) is in agreement 
with experiments for both processes. The frequency of 



cessations can be calculated based on an analogy to the 
Arrhenius-Kramers problem. The azimuthal dynamics 
during cessations can be understood in terms of the an- 
gular momentum of the LSC which suppresses azimuthal 
motion driven by turbulent fluctuations. The model also 
can reproduce other results for the LSC dynamics in good 
agreement with experiments; these are to be studied thor- 
oughly in future work. 
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